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Abstract 

The  funding  from  this  grant  was  utilized  to  further  develop  a  new  approach  for  the  elec¬ 
tronic  structure  of  materials.  The  method  formulates  the  Kohn-Sham  equations  of  Density 
Functional  Theory  directly  in  real  space  with  a  high  order  Finite  Difference  approach.  The 
resulting  equations  were  solved  using  the  linear  scaling  multigrid  algorithm  developed  by 
Brandt  and  coworkers.  Multigrid  techniques  were  used  to  solve  both  the  self  consistent 
eigenvalue  equations  and  the  Poisson  equation  for  the  electrostatic  potential  at  each  step 
of  iterations.  Accurate  numerical  results  were  obtained  for  finite  and  periodic  electrostatic 
problems  and  for  the  eigenvalue  equations  for  many  electron  atoms  and  simple  molecules. 
Recently,  conservative  grid  equations  have  been  developed  so  that  grid  refinement  strategies 
can  be  employed.  This  allows  one  to  perform  extensive  numerical  work  selectively  in  regions 
of  high  electron  density.  The  new  method  should  have  wide  applications  for  numerical  stud¬ 
ies  of  complex  and  disordered  materials  which  require  a  quantum  mechanical  treatment  for 
many  atoms. 


Final  Report 

Many  computations  in  chemistry  require  treatments  which  go  beyond  effective  pair  or 
triplet  intermolecular  potentials.  Examples  are  chemical  reactions  in  solution,  plasma  chem- 
istr\%  glasses  with  complex  covalent  networks,  high  pressure  phase  diagrams,  interactions 
of  molecules  with  metal  surfaces,  responses  of  large  systems  to  electromagnetic  fields,  etc. 
Quantum  mechanical  computations  on  large  systems  have  been  limited  by  the  scaling  proper¬ 
ties  of  existing  algorithms.  In  recent  years  extensive  effort  has  gone  into  developing  methods 
which  result  in  better  scaling  properties.  With  the  support  from  AFOSR,  we  have  developed 
a  new  approach  which  utilizes  a  multigrid  (MG)  solver^’^  for  both  the  Kohn-Sham  equations 
and  the  Poisson  equation.^”®  The  method  scales  in  computer  time  linearly  with  the  number 
of  electrons  if  the  orbitals  can  be  localized  in  real  space  (and  still  yield  a  realistic  descrip¬ 
tion  of  the  ground  state  electron  density).  Our  calculations  so  far  have  been  all  electron 
calculations  with  bare  nuclear  charges. 

In  our  first  study, we  developed  a  multiscale  approach  for  solving  the  Kohn-Sham  equa¬ 
tions  in  real  space.  A  high  order  (8th)  finite  difference  (FD)  representation  was  used  for  the 
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kinetic  energy  operator  in  order  to  obtain  numerical  accuracy  with  reasonable  sized  grids. 
The  nucleus  was  represented  as  a  discrete  charge  on  the  grid  with  spacial  extent  that  of  the 
finest  scale  grid  size  h.  The  Kohn-Sham  equations  were  solved  first  on  the  coarse  scale;  then 
the  solution  was  interpolated  to  the  next  finer  scale  and  iterations  were  initiated  there.  This 
method  does  not  guarantee  a  linear  scaling  solution  due  to. the  phenomenon  of  aliasing  of 
high  frequency  modes  on  the  coarse  scale.  However,  we  obtained  accurate  numerical  results 
in  all  electron  calculations  with  uniform  grids,  and  the  multiscale  method  led  to  substantial 
accelerations  toward  the  solution.  This  reflects  the  presence  of  critical  slowing  down  (CSD) 
in  the  solution  process  and  indicates  the  importance  of  multiple  length  scales  in  the  elec¬ 
tronic  structure  problem.  By  just  iterating  on  the  fine  scale  alone,  it  appears  one  is  not  even 
guaranteed  of  approaching  the  correct  solution  in  a  reasonable  amount  of  time. 

In  a  second  paper,^  we  presented  our  high  order  MG  method  for  solving  the  Poisson 
equation  which  is  required  at  each  self  consistency  iteration  in  the  Kohn-Sham  equations. 
We  tested  the  method  on  several  finite  and  periodic  problems  to  determine  its  efficiency 
and  accuracy.  We  obtained  accurate  results  for  a  simple  model  for  atoms  with  a  discrete 
nuclear  charge  at  the  origin  and  a  balancing  electron  distribution.  Accurate  results  were 
also  obtained  for  periodic  lattice  computations  in  3-d.  This  means  the  method  can  handle 
long  ranged  Coulomb  problems  without  requiring  implementation  of  the  Ewald  method.  By 
placing  the  nuclei  on  the  grid  as  discrete  charges,  the  electrons  and  nuclei  are  on  an  equal 
footing  in  terms  of  solution  of  the  Poisson  equation.  The  linear  scaling  attributes  of  the 
algorithm  were  demonstrated  numerically.  The  solver  was  used  to  generate  an  accurate 
electrostatic  map  of  a  model  zeolite  solid.  In  a  subsequent  report,®  similar  computations 
were  performed  for  a  4096  ion  solid.  The  complete  solution  of  the  Poisson  equation  on  a 
128®  grid  starting  from  <;/)(r)=0  everywhere  took  3  seconds  on  the  OSC  Cray-YMP.  We  are 
currently  using  the  high  order  Poisson  solvei'  to  carry  out  Monte  Carlo  simulations  of  ionic 
fluids  to  test  its  viability  as  a  general  method  for  electrostatic  computations  in  molecular 
simulations.  ..  . 

Recently,  we  have  successfully  incorporated  the  nonlinear  full  approximation  scheme 
(F.4S)  MG  eigenvalue  solver  of  Brandt,  et  al^  into  our  method.^  This  work  was  presented 
at  the  International  Paris  DFT  Meeting  of  1995.  With  this  algorithm,  we  obtained  accurate 
ground  state  results  for  the  Ne  atom  (all  electrons)  with  on  the  order  of  only  10  fine  scale 
iterations.  The  full  MG  solver  outperforms  the  nested  procedure,  although  the  nested  itera¬ 
tion  does  yield  comparably  accurate  results.  The  solution  of  the  Kohn-Sham  equations  is  a 
nonlinear  problem,  which  requires  an  MG  method  going  bej^ond  the  standard  approach  typ¬ 
ical  for  linear  elliptic  problems  such  as  the  Poisson  equation.  We  summarized  our  approach, 
numerical  results,  limitations,  and  future  plans  in  this  paper. 

In  recent  months  we  have  begin  to  examine  grid  refinement  strategies,  which  are  required 
to  obtain  accurate  numerical  results  with  reasonable  grid  sizes. We  have  utlized  the  tech¬ 
niques  of  Bai  and  Brandt^^  for  solving  the  Poisson  equation  in  regions  where  the  potential 
varies  rapidly,  namely  near  the  nucleus.  Their  formulation  was  developed  for  2nd  order 
equations,  and  we  are  currently  working  in  collaboration  with  Achi  Brandt  on  a  higher  order 
method.  Conservative  equations  are  required  because,  when  the  difference  equations  are  set 
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up  for  the  coarse  scale  iterations,  the  ‘defect  correction’  r  leads  to  a  nonconserved  source 
within  the  finite  refinement  zone.  This  can  pollute  the  solution  over  long  ranges  for  singular 
sources.  We  have  implemented  for  the  first  time  the  Bai  and  Brandt  technique  for  three 
dimensional  finite  and  periodic  domains  and  singular  sources,  and  have  obtained  promising 
numerical  results  with  their  2nd  order  method.  The  higher  order  formulation  will  be  crucial 
in  solving  real  world  problems  on  grids  of  approachable  sizes.  These  techniques  should  have 
applications  in  several  areas  of  computational  physics  and  chemistry.  Our  focus  wdll  be  to 
use  the  refinement  methods  in  the  region  around  the  atomic  nucleus  to  obtain  increased 
resolution  in  the  core. 

In  a  related  project,  we  have  developed  a  novel  approach  for  computation  of  chemi¬ 
cal  potentials  and  phase  equlibria  in  dense  fluids. The  method  uses  the  Kirkwood-Buff 
theory  of  solutions  and  an  inversion  procedure  developed  by  Ben-Naim  to  compute  ther¬ 
modynamic  derivative  information.  This  information  is  used  to  extrapolate  the  chemical 
potentials  around  a  given  state  point.  The  phase  diagram  of  high  pressure  He/H2  fluid  mix¬ 
tures  .was  successfully  computed.  Our  intention  is  to  go  to  even  higher  pressures  (near  1 
I\fegabar)  where  the  hydrogen  component  can  become  metallic.  The  phase  diagram  of  the 
mixture  is  unknown  there,  and  has  an  impact  on  our  understanding  of  the  evolution  of  the 
planet  Jupiter.  The  multigrid  electronic  structure  solver  will  be  used  in  conjunction  with 
the  statistical  mechanical  methods  to  attempt  for  the  first  time  to  compute  a  physical  phase 
diagram  from  ab  initio  quantum  mechanics.  If  the  work  is  successful  it  will  pave  the  way  for 
simulations  of  more  chemically  complex  mixtures,  such  as  occur  in  the  interior  of  Earth. 

The  funding  from  this  AFOSR  grant  was  used  to  partially  support  two  fine  postdoctoral 
researchers,  Dr.  Karthik  l3'er  and  Dr.  .A.nping  Liu,  for  several  months  while  they  worked  on 
the  FAS  eigenvalue  problem  and  related  statistical  mechanical  computations.  The  research 
is  currently  funded  by  an  NSF  grant  and  further  support  will  be  sought  from  AFOSR  in  the 
near  future.  Reprints  of  our  research  papers  our  enclosed. 
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